rm(list = ls())

#setwd("~/Dropbox/Energiewende_Project/All_Replication_Files")
setwd("C:/Users/chris/Dropbox/Energiewende_Project/All_Replication_Files")

#Loading in Packages ----------------------------------------------------------------------------------

library(dplyr)
library(sandwich)
library(lmtest)
library(margins)
library(ggplot2)
library(scales)
library(gridExtra)
library(readr)
library(stargazer)
library(reshape2) 
library(stringr)
library(cregg)
library(kableExtra)
library(knitr)
library(modelsummary)

####################################################################################################################################
############################################Coal Exit Experiment####################################################################
####################################################################################################################################
###########Prep Data:####
data <- read_csv("data.coal.cleaned.final.csv", col_types = cols(.default = "?"))
weights <- read_csv("weights_data.csv")

data <- merge(data, weights, by=c("tic"), all.x=T)

data$coalexit <- relevel(as.factor(data$coalexit), "Control")

###########Helper Functions:#####
#main regression:
reg_coal <- function(data, subgroup){
  
  datac <- data[data[,subgroup] == 1,]
  fmla <- function(x){
    as.formula(paste0(x, "~ coalexit"))
  }
  mods.coal = lapply(y.variables, fmla)
  
  mod.fit = function(fmla){
    lm(fmla, data=datac, weights=datac$weights)
  }
  
  fits.coal = lapply(mods.coal, mod.fit)

  return(fits.coal)
}

#with covariates

reg_coal_cov <- function(data, subgroup){
  
  datac <- data[data[,subgroup] == 1,]
  
  fmla <- function(x){
    as.formula(paste0(x, "~ coalexit + region + age + gender + incomequintiles "))
  }
  mods.coal = lapply(y.variables, fmla)
  
  mod.fit = function(fmla){
    lm(fmla, data=datac, weights=datac$weights)
  }
  
  fits.coal = lapply(mods.coal, mod.fit)
  
  return(fits.coal)
}

#with heterogeneity interaction terms:


reg_coal_int<- function(data,inter, subgroup){
  
  datac <- data[data[,subgroup] == 1,]

    fmla <- function(x){
      as.formula(paste0(x, "~ coalexit +", inter, "+ coalexit*", inter))
    }
    mods.coal = lapply(y.variables, fmla)

    mod.fit = function(fmla){
      lm(fmla, data=datac, weights=datac$weights)
    }

  fits.coal = lapply(mods.coal, mod.fit)
   
  return(fits.coal)
}

reg_coal_int_cov<- function(data,inter, subgroup){
  
  datac <- data[data[,subgroup] == 1,]

  
  fmla <- function(x){
    as.formula(paste0(x, "~ coalexit + region + age + gender + incomequintiles+", inter, "+ coalexit*", inter))
  }
  mods.coal = lapply(y.variables, fmla)
  
  mod.fit = function(fmla){
    lm(fmla, data=datac, weights=datac$weights)
  }
  
  fits.coal = lapply(mods.coal, mod.fit)
  

  return(fits.coal)
}


##robustness to other attitudes


reg_coal_int_rob <- function(data,inter, subgroup){
  
  datac <- data[data[,subgroup] == 1,]
  
  fmla <- function(x){
    as.formula(paste0(x, "~ coalexit + region + age + gender + incomequintiles+", inter, "+ coalexit*", inter, "+ climatehumans*coalexit+climate_policyn*coalexit+ left_right_4*coalexit+ party_2021_2*coalexit"))
  }
  mods.coal = lapply(y.variables, fmla)
  
  mod.fit = function(fmla){
    lm(fmla, data=datac, weights=datac$weights)
  }
  
  fits.coal = lapply(mods.coal, mod.fit)
  
  return(fits.coal)
}


#data prep for plotting:
data.prep <- function(reg){
  regsum <- summary(reg)
  coef <- as.matrix(regsum$coefficients[2:5,1:2])  
  coef <- as.data.frame(cbind(c("Compensation for Consumers Treatment", "Compensation for Investors Treatment", 
                                "Compensation for Regions Treatment", "Compensation for Workers Treatment"), coef))
  coef$Estimate <- as.numeric(paste(coef$Estimate))
  coef$`Std. Error`<- as.numeric(paste(coef$`Std. Error`))
  coef$confupper <- coef$Estimate + coef$`Std. Error`*1.96
  coef$conflower <- coef$Estimate - coef$`Std. Error`*1.96
  coef$V1 <- as.factor(coef$V1)
  coef$V1 <- factor(coef$V1,
                    levels = levels(coef$V1)[c(1,2,3,4,5)])
  return(coef)
}

#regression tables:

#plot results:
plot_basic <- function(regout, fileadd, dvadd){
  data.temp <- data.prep(regout)
  
  pl <-  ggplot(data.temp)+
    geom_point(aes(x = Estimate, y = V1))+
    geom_segment(aes(x = conflower, xend = confupper,
                     y = V1, yend = V1))+
    geom_vline(xintercept = 0, lty = "dotted")+
    theme_bw()+
    ggtitle(paste0("Treatment Effects Coal Vignette, \n", dvadd, ", \n Weighted Regression"))
  
  
  pl <- pl + theme(text = element_text(size=16),
                   plot.title = element_text(size=18, hjust = 0.5),
                   axis.text=element_text(size=16),
                   axis.title = element_text(size=16))+
    scale_y_discrete(labels = function(x) str_wrap(x, width = 15))+
    xlab("Average Treatment Effects")+
    ylab("")
  
  print(pl)
  
  ggsave(filename = paste0("Graphs/coal_vignette.weighted.", fileadd, ".png"),
         plot = pl,
         width = 12, height = 8)
}


#heterogeneity analysis:
drawplots_coal <- function(res, n.res, colors, legendtitle ,addtitle=NULL, fileadd = NULL){
  for(k in 1:4){
    plot.datac <- lapply(res,"[[",k)
    
    plot.datap <- lapply(plot.datac, data.prep)
    
    data.temp <- as.data.frame(plot.datap[[1]])
    
    pl <-  ggplot(data.temp)+
      geom_point(aes(x = Estimate, y = V1, color = paste0(colors[1] , ", n= ", n.res[[1]]), shape = paste0(colors[1] , ", n= ", n.res[[1]])))+
      geom_segment(aes(x = conflower, xend = confupper,
                       y = V1, yend = V1, color = paste0(colors[1] , ", n= ", n.res[[1]])))+
      geom_vline(xintercept = 0, lty = "dotted")+
      theme_bw()+
      ggtitle(paste0("Treatment Effects Coal Vignette", addtitle, "\n ", dependent.variables[k], ", \n Weighted Regression"))
    
    
    
    
    if(length(res) > 1){
      data.temp <- as.data.frame(plot.datap[[2]])
      pl <- pl +  geom_point(data = data.temp, aes(x = Estimate, y = V1, color = paste0(colors[2] , ", n= ", n.res[[2]]), shape = paste0(colors[2] , ", n= ", n.res[[2]])), position = position_nudge(y = -0.1))+
        geom_segment(data = data.temp,aes(x = conflower, xend = confupper,
                                          y = V1, yend = V1, color = paste0(colors[2] , ", n= ", n.res[[2]])), position = position_nudge(y = -0.1))
    }
    
    if(length(res) > 2){
      data.temp <- as.data.frame(plot.datap[[3]])
      pl <- pl +  geom_point(data = data.temp, aes(x = Estimate, y = V1, color = paste0(colors[3] , ", n= ", n.res[[3]])), position = position_nudge(y = +0.1))+
        geom_segment(data = data.temp,aes(x = conflower, xend = confupper,
                                          y = V1, yend = V1, color = paste0(colors[3] , ", n= ", n.res[[3]])), position = position_nudge(y = +0.1))
    }
    
    if(length(res) > 3){
      data.temp <- as.data.frame(plot.datap[[4]])
      pl <- pl +  geom_point(data = data.temp, aes(x = Estimate, y = V1, color = paste0(colors[4] , ", n= ", n.res[[4]])), position = position_nudge(y = -0.2))+
        geom_segment(data = data.temp,aes(x = conflower, xend = confupper,
                                          y = V1, yend = V1, color = paste0(colors[4] , ", n= ", n.res[[4]])), position = position_nudge(y = -0.2))
    }
    
    
    if(length(res) > 4){
      data.temp <- as.data.frame(plot.datap[[5]])
      pl <- pl +  geom_point(data = data.temp, aes(x = Estimate, y = V1, color = paste0(colors[5] , ", n= ", n.res[[5]])), position = position_nudge(y = +0.2))+
        geom_segment(data = data.temp,aes(x = conflower, xend = confupper,
                                          y = V1, yend = V1, color = paste0(colors[5] , ", n= ", n.res[[5]])), position = position_nudge(y = +0.2))
    }
    
    if(length(res) > 5){
      data.temp <- as.data.frame(plot.datap[[6]])
      pl <- pl +  geom_point(data = data.temp, aes(x = Estimate, y = V1, color = paste0(colors[6] , ", n= ", n.res[[6]])), position = position_nudge(y = -0.3))+
        geom_segment(data = data.temp,aes(x = conflower, xend = confupper,
                                          y = V1, yend = V1, color = paste0(colors[6] , ", n= ", n.res[[6]])), position = position_nudge(y = -0.3))
    }
    
    pl <- pl + labs(color='Group') + theme(text = element_text(size=16),
                                           plot.title = element_text(size=18, hjust = 0.5),
                                           axis.text=element_text(size=16),
                                           axis.title = element_text(size=16),
                                           legend.position = "bottom",
                                           legend.text = element_text(size=16))+
      scale_y_discrete(labels = function(x) str_wrap(x, width = 15))+
      xlab("Average Treatment Effects")+
      ylab("") + 
      guides(color=guide_legend(title=legendtitle), 
             shape = guide_legend(title=legendtitle))+
      scale_color_grey(start = 0.1, end = 0.5)
    
    print(pl)
    
    ggsave(filename = paste("Graphs/coal_vignette.weighted.", fileadd, file.names.y[k],".png", sep=""),
           plot = pl,
           width = 12, height = 8)
    
  }
  
}


###########Define Terms and Vectors:#####
y.variables <- c("coalexitpremax", "coalexitsooner", "E3_exitsupport","coalexitsupport")

dependent.variables <- c("2050 - Preferred Year of Exit", "Coal Exit Sooner than 2038", "Support for Full Exit, Continuous", "Support for Full Exit, Binary")

file.names.y <- c("maxyear", "sooner", "supportc", "supportb")

supportinteractions <- c("highsupportinterventionecon", "highsupportinterventionfrail", "highsupportinterventionworkers")


###########Basic Regression:####

fits.coal <- reg_coal_cov(data, "all")


#drawing plots of four main regressions:

for(i in 1:4){
  plot_basic(fits.coal[[i]], file.names.y[i], dependent.variables[i])
}


####heterogeneity plots with split samples

plot.data <- list()
n.subsets <- list()
subsets <-  c("highsupportintervention", "lowsupportintervention")

for(i in 1:length(subsets)){
  regressionsout <- reg_coal_cov(data, subsets[i])
  plot.data[[i]] <- regressionsout
  n.subsets[[i]] <- length(regressionsout[[1]]$fitted.values)
}

drawplots_coal(plot.data, n.subsets, c("High", "Low") , paste0("Belief State Intervention"),  " by Belief in State Intervention", "govtint") 


##########Creating Regression Tables for Paper:###############

#table 1: baseline and interaction with belief in state management

for(i in 1:4){
  reg.1 <- reg_coal(data, "all")[[i]]
  
  reg.2 <-  reg_coal_cov(data, "all")[[i]]
  
  reg.3 <- reg_coal_int(data, "highsupportintervention", "all")[[i]]
  
  reg.4 <- reg_coal_int_cov(data, "highsupportintervention", "all")[[i]]
  
  reg.5 <- reg_coal_int(data, "supportintervention", "all")[[i]]
  
  reg.6 <- reg_coal_int_cov(data, "supportintervention", "all")[[i]]
  
  
  stargazer(reg.2,reg.4,  reg.6,
            keep=c("coalexitCompensation for Consumers",
                   "coalexitCompensation for Investors",
                   "coalexitCompensation for Regions",
                   "coalexitCompensation for Workers",
                   "highsupportintervention",
                   "coalexitCompensation for Consumers:highsupportintervention",
                   "coalexitCompensation for Investors:highsupportintervention",
                   "coalexitCompensation for Regions:highsupportintervention",
                   "coalexitCompensation for Workers:highsupportintervention",
                   "supportintervention",
                   "coalexitCompensation for Consumers:supportintervention",
                   "coalexitCompensation for Investors:supportintervention",
                   "coalexitCompensation for Regions:supportintervention",
                   "coalexitCompensation for Workers:supportintervention"),
            covariate.labels =c("Compensation for Consumers","Compensation for Investors",
                                "Compensation for Regions","Compensation for Workers",
                                "Belief State Int. Bin.", "Comp. Consumers * Belief State Int. Bin.",
                                "Comp. Investors * Belief State Int. Bin.",
                                "Comp. Regions * Belief State Int. Bin.",
                                "Comp. Workers * Belief State Int. Bin.",
                                "Belief State Int. Cont.", "Comp. Consumers * Belief State Int. Cont.",
                                "Comp. Investors * Belief State Int. Cont.",
                                "Comp. Regions * Belief State Int. Cont.",
                                "Comp. Workers * Belief State Int. Cont."),
            add.lines = list(c("Demographic Controls", "Yes",  "Yes",  "Yes")),
            font.size="tiny",
            dep.var.labels = dependent.variables[i],
            out=paste0("Tables/coal_vignette.weighted.",file.names.y[i],".tex")
  )
  
  
}
# 
# ###full table with demographics
# 
# for(i in 1:4){
#   reg.1 <- reg_coal(data, "all")[[i]]
#   
#   reg.2 <-  reg_coal_cov(data, "all")[[i]]
#   
#   reg.3 <- reg_coal_int(data, "highsupportintervention", "all")[[i]]
#   
#   reg.4 <- reg_coal_int_cov(data, "highsupportintervention", "all")[[i]]
#   
#   reg.5 <- reg_coal_int(data, "supportintervention", "all")[[i]]
#   
#   reg.6 <- reg_coal_int(data, "supportintervention", "all")[[i]]
#   
#   
#   table.1.full <- modelsummary(list(reg.1, reg.2, reg.3, reg.4, reg.5, reg.6),
#                                stars=T,
#                                coef_rename = c("coalexitCompensation for Consumers Treatment" = "Compensation for Consumers",
#                                                "coalexitCompensation for Investors Treatment" = "Compensation for Investors",
#                                                "coalexitCompensation for Regions" = "Compensation for Regions",
#                                                "coalexitCompensation for Workers" = "Compensation for Workers",
#                                                "highsupportintervention"= "Belief St. Int. Bin.",
#                                                "supportintervention" = "Belief in St. Int. Cont."
#                                ),
#                                output="latex")
#   
#   header.t1 <- paste("DV:", dependent.variables[i])
#   
#   tableheader.t1 <- c(c(" " = 1,  header.t1 = 6))
#   
#   # set vector names 
#   names(tableheader.t1) <- c(" ", header.t1)
#   
#   
#   regtable.t1 <- table.1.full %>% kableExtra::add_header_above(header=tableheader.t1) 
#   
#   
#   save_kable(regtable.t1, paste0("Tables/coal_vignette_full_weighted.",file.names.y[i],".tex"))
#   
#   
# }
# 
# 
# ##table 1 repeated for only those who passed the mock vignette
# 
# for(i in 1:4){
#   reg.1.mv <- reg_coal(data, "coalattention")[[i]]
#   
#   reg.2.mv <-  reg_coal_cov(data, "coalattention")[[i]]
#   
#   reg.3.mv <- reg_coal_int(data, "highsupportintervention", "coalattention")[[i]]
#   
#   reg.4.mv <- reg_coal_int_cov(data, "highsupportintervention", "coalattention")[[i]]
#   
#   reg.5.mv <- reg_coal_int(data, "supportintervention", "coalattention")[[i]]
#   
#   reg.6.mv <- reg_coal_int(data, "supportintervention", "coalattention")[[i]]
#   
#   
#   stargazer(reg.1.mv, reg.2.mv, reg.3.mv, reg.4.mv, reg.5.mv, reg.6.mv,
#             keep=c("coalexitCompensation for Consumers",
#                    "coalexitCompensation for Investors",
#                    "coalexitCompensation for Regions",
#                    "coalexitCompensation for Workers",
#                    "highsupportintervention",
#                    "coalexitCompensation for Consumers:highsupportintervention",
#                    "coalexitCompensation for Investors:highsupportintervention",
#                    "coalexitCompensation for Regions:highsupportintervention",
#                    "coalexitCompensation for Workers:highsupportintervention",
#                    "supportintervention",
#                    "coalexitCompensation for Consumers:supportintervention",
#                    "coalexitCompensation for Investors:supportintervention",
#                    "coalexitCompensation for Regions:supportintervention",
#                    "coalexitCompensation for Workers:supportintervention"),
#             covariate.labels =c("Compensation for Consumers","Compensation for Investors",
#                                 "Compensation for Regions","Compensation for Workers",
#                                 "Belief State Int. Bin.", "Comp. Consumers * Belief State Int. Bin.",
#                                 "Comp. Investors * Belief State Int. Bin.",
#                                 "Comp. Regions * Belief State Int. Bin.",
#                                 "Comp. Workers * Belief State Int. Bin.",
#                                 "Belief State Int. Cont.", "Comp. Consumers * Belief State Int. Cont.",
#                                 "Comp. Investors * Belief State Int. Cont.",
#                                 "Comp. Regions * Belief State Int. Cont.",
#                                 "Comp. Workers * Belief State Int. Cont."),
#             add.lines = list(c("Demographic Controls", "No" , "Yes", "No" , "Yes", "No" , "Yes")),
#             font.size="tiny",
#             dep.var.labels = dependent.variables[i],
#             out=paste0("Tables/coal_vignette_mv_passed.weighted.",file.names.y[i],".tex")
#   )
#   
#   
# }
# 
# ###full table with demographics
# 
# for(i in 1:4){
#   reg.1.mv <- reg_coal(data, "coalattention")[[i]]
#   
#   reg.2.mv <-  reg_coal_cov(data, "coalattention")[[i]]
#   
#   reg.3.mv <- reg_coal_int(data, "highsupportintervention", "coalattention")[[i]]
#   
#   reg.4.mv <- reg_coal_int_cov(data, "highsupportintervention", "coalattention")[[i]]
#   
#   reg.5.mv <- reg_coal_int(data, "supportintervention", "coalattention")[[i]]
#   
#   reg.6.mv <- reg_coal_int(data, "supportintervention", "coalattention")[[i]]
#   
#   
#   table.1.full <- modelsummary(list(reg.1, reg.2, reg.3, reg.4, reg.5, reg.6),
#                                stars=T,
#                                coef_rename = c("coalexitCompensation for Consumers Treatment" = "Compensation for Consumers",
#                                                "coalexitCompensation for Investors Treatment" = "Compensation for Investors",
#                                                "coalexitCompensation for Regions" = "Compensation for Regions",
#                                                "coalexitCompensation for Workers" = "Compensation for Workers",
#                                                "highsupportintervention"= "Belief St. Int. Bin.",
#                                                "supportintervention" = "Belief in St. Int. Cont."
#                                ),
#                                output="latex")
#   
#   header.t1 <- paste("DV:", dependent.variables[i])
#   
#   tableheader.t1 <- c(c(" " = 1,  header.t1 = 6))
#   
#   # set vector names 
#   names(tableheader.t1) <- c(" ", header.t1)
#   
#   
#   regtable.t1 <- table.1.full %>% kableExtra::add_header_above(header=tableheader.t1) 
#   
#   
#   save_kable(regtable.t1, paste0("Tables/coal_vignette_attention_full_weighted.",file.names.y[i],".tex"))
#   
#   
# }

# #table 2, interactions with different subsets of government policy:
# 
# list.table2 <- list()
# 
# 
# for(i in 1:4){
#   
#   for(inter in supportinteractions){
#     list.table2[[inter]] <- reg_coal_int(data, inter, "all")[[i]]
#   }
#   
#   tab1 <- list.table2[["highsupportinterventionecon"]]
#   tab2 <- list.table2[["highsupportinterventionfrail"]]
#   tab3 <- list.table2[["highsupportinterventionworkers"]]
#   
#   stargazer(tab1, tab2, tab3,
#             keep=c("coalexitCompensation for Consumers",
#                    "coalexitCompensation for Investors",
#                    "coalexitCompensation for Regions",
#                    "coalexitCompensation for Workers",
#                    "highsupportinterventionecon",
#                    "coalexitCompensation for Consumers:highsupportinterventionecon",
#                    "coalexitCompensation for Investors:highsupportinterventionecon",
#                    "coalexitCompensation for Regions:highsupportinterventionecon",
#                    "coalexitCompensation for Workers:highsupportinterventionecon",
#                    "highsupportinterventionfrail",
#                    "coalexitCompensation for Consumers:highsupportinterventionfrail",
#                    "coalexitCompensation for Investors:highsupportinterventionfrail",
#                    "coalexitCompensation for Regions:highsupportinterventionfrail",
#                    "coalexitCompensation for Workers:highsupportinterventionfrail",
#                    "highsupportinterventionworkers",
#                    "coalexitCompensation for Consumers:highsupportinterventionworkers",
#                    "coalexitCompensation for Investors:highsupportinterventionworkers",
#                    "coalexitCompensation for Regions:highsupportinterventionworkers",
#                    "coalexitCompensation for Workers:highsupportinterventionworkers"
#             ),
#             covariate.labels =c("Compensation for Consumers","Compensation for Investors",
#                                 "Compensation for Regions","Compensation for Workers",
#                                 "Belief State Int. Economy Bin.", "Comp. Consumers * Belief State Int. Economy Bin.",
#                                 "Comp. Investors * Belief State Int. Economy Bin.",
#                                 "Comp. Regions * Belief State Int. Economy Bin.",
#                                 "Comp. Workers * Belief State Int. Economy Bin.",
#                                 "Belief State Int. Frail Bin.", "Comp. Consumers * Belief State Int. Frail Bin.",
#                                 "Comp. Investors * Belief State Int. Frail Bin.",
#                                 "Comp. Regions * Belief State Int. Frail Bin.",
#                                 "Comp. Workers * Belief State Int. Frail Bin.",
#                                 "Belief State Int. Workers Bin.", "Comp. Consumers * Belief State Int. Workers Bin.",
#                                 "Comp. Investors * Belief State Int. Workers Bin.",
#                                 "Comp. Regions * Belief State Int. Workers Bin.",
#                                 "Comp. Workers * Belief State Int. Workers Bin."),
#             font.size="tiny",
#             dep.var.labels = dependent.variables[i],
#             out=paste0("Tables/coal_vignette_sub_policy.weighted.",file.names.y[i],".tex")
#   )
#   
#   
# }
# 
# #only those who passed the mock vignette before coal experiment
# 
# list.table2 <- list()
# 
# for(i in 1:4){
#   
#   for(inter in supportinteractions){
#     list.table2[[inter]] <- reg_coal_int(data, inter, "coalattention")[[i]]
#   }
#   
#   tab1 <- list.table2[["highsupportinterventionecon"]]
#   tab2 <- list.table2[["highsupportinterventionfrail"]]
#   tab3 <- list.table2[["highsupportinterventionworkers"]]
#   
#   stargazer(tab1, tab2, tab3,
#             keep=c("coalexitCompensation for Consumers",
#                    "coalexitCompensation for Investors",
#                    "coalexitCompensation for Regions",
#                    "coalexitCompensation for Workers",
#                    "highsupportinterventionecon",
#                    "coalexitCompensation for Consumers:highsupportinterventionecon",
#                    "coalexitCompensation for Investors:highsupportinterventionecon",
#                    "coalexitCompensation for Regions:highsupportinterventionecon",
#                    "coalexitCompensation for Workers:highsupportinterventionecon",
#                    "highsupportinterventionfrail",
#                    "coalexitCompensation for Consumers:highsupportinterventionfrail",
#                    "coalexitCompensation for Investors:highsupportinterventionfrail",
#                    "coalexitCompensation for Regions:highsupportinterventionfrail",
#                    "coalexitCompensation for Workers:highsupportinterventionfrail",
#                    "highsupportinterventionworkers",
#                    "coalexitCompensation for Consumers:highsupportinterventionworkers",
#                    "coalexitCompensation for Investors:highsupportinterventionworkers",
#                    "coalexitCompensation for Regions:highsupportinterventionworkers",
#                    "coalexitCompensation for Workers:highsupportinterventionworkers"
#             ),
#             covariate.labels =c("Compensation for Consumers","Compensation for Investors",
#                                 "Compensation for Regions","Compensation for Workers",
#                                 "Belief State Int. Economy Bin.", "Comp. Consumers * Belief State Int. Economy Bin.",
#                                 "Comp. Investors * Belief State Int. Economy Bin.",
#                                 "Comp. Regions * Belief State Int. Economy Bin.",
#                                 "Comp. Workers * Belief State Int. Economy Bin.",
#                                 "Belief State Int. Frail Bin.", "Comp. Consumers * Belief State Int. Frail Bin.",
#                                 "Comp. Investors * Belief State Int. Frail Bin.",
#                                 "Comp. Regions * Belief State Int. Frail Bin.",
#                                 "Comp. Workers * Belief State Int. Frail Bin.",
#                                 "Belief State Int. Workers Bin.", "Comp. Consumers * Belief State Int. Workers Bin.",
#                                 "Comp. Investors * Belief State Int. Workers Bin.",
#                                 "Comp. Regions * Belief State Int. Workers Bin.",
#                                 "Comp. Workers * Belief State Int. Workers Bin."),
#             font.size="tiny",
#             dep.var.labels = dependent.variables[i],
#             out=paste0("Tables/coal_vignette_sub_policy_mv.weighted.",file.names.y[i],".tex")
#   )
#   
#   
# }

#split sample tables

list.tablesplit.h <- reg_coal_cov(data, "highsupportintervention")
list.tablesplit.l <- reg_coal_cov(data, "lowsupportintervention")

list.tablesplit.h.a <- reg_coal_cov(data[data$coalattention==1,], "highsupportintervention")
list.tablesplit.l.a <- reg_coal_cov(data[data$coalattention==1,], "lowsupportintervention")

for(i in 1:4){
  
  table_split <- modelsummary(list(list.tablesplit.h[[i]], list.tablesplit.l[[i]], list.tablesplit.h.a[[i]], list.tablesplit.l.a[[i]]),
                              stars=T,
                              coef_rename = c("coalexitCompensation for Consumers" = "Compensation for Consumers",
                                              "coalexitCompensation for Investors" = "Compensation for Investors",
                                              "coalexitCompensation for Regions" = "Compensation for Regions",
                                              "coalexitCompensation for Workers" = "Compensation for Workers"),
                              coef_omit="age|region|gender|income",
                              output="latex")
  
  header.t1 <- paste("DV:", dependent.variables[i])
  
  tableheader.ts1 <- c(c(" " = 1,  header.t1 = 4 ))
  
  # set vector names 
  names(tableheader.ts1) <- c(" ", header.t1)
  
  high_a <- "Passed Attention"
  all <- "All"
  
  tableheader.ts2 <- c(" " = 1, all=2, high_a = 2)
  names(tableheader.ts2) <- c(" ", all, high_a)
  
  high_b <- "High Belief St. Int."
  low_b <- "Low Belief St. Int."
  
  tableheader.ts3 <- c(" " = 1, high_b =1, low_b = 1, high_b=1, low_b=1)
  names(tableheader.ts3) <- c(" ", high_b, low_b, high_b, low_b)
  
  regtable.ts1 <- table_split %>% kableExtra::add_header_above(header=tableheader.ts3)%>% kableExtra::add_header_above(header=tableheader.ts2) %>% kableExtra::add_header_above(header=tableheader.ts1) 
  
  
  save_kable(regtable.ts1, paste0("Tables/coal_vignette_split_weighted.",file.names.y[i],".tex"))
  
}



###Table 3: Robust interactions with interactions with other attitude variables

for(i in 1:4){
  
  
  header.t1 <- paste("DV:", dependent.variables[i])
  
  tableheader.tr1 <- c(c(" " = 1,  header.t1 = 1))
  
  reg.1.r <- reg_coal_int_rob(data, "highsupportintervention", "all")[[i]]
  
  
  
  table.rf1 <- modelsummary(list(reg.1.r), stars=T,
                            gof_omit="AIC|BIC",
                            coef_rename = c("coalexitCompensation for Consumers" = "Compensation for Consumers",
                                            "coalexitCompensation for Investors" = "Compensation for Investors",
                                            "coalexitCompensation for Regions" = "Compensation for Regions",
                                            "coalexitCompensation for Workers" = "Compensation for Workers",
                                            "party_2021_2" = "Party Vote 2021",
                                            "left_right_4" = "Left-Right Leaning",
                                            "climate_policyn"= "Climate Policy Support",
                                            "climatehumans" = "Belief in MM CC",
                                            "highintervention" = "Belief in St. Int. Bin.",
                                            "supportintervention" = "Belief in St. Int. Cont.",
                                            "region" = "Region "),
                            estimate  = "{estimate} ({std.error})",
                            statistic = NULL,
                            output="latex")
  
  
  
  regtable.rf1 <- table.rf1 %>% kableExtra::add_header_above(header=tableheader.tr1) 
  
  
  save_kable(regtable.rf1, paste0("Tables/coal_vignette_robust1.weighted.",file.names.y[i],".tex"))
  
  
  
  reg.2.r <- reg_coal_int_rob(data, "supportintervention", "all")[[i]]
  
  table.rf2 <- modelsummary(list(reg.2.r), stars=T,
                            gof_omit="AIC|BIC",
                            coef_rename = c("coalexitCompensation for Consumers Treatment" = "Compensation for Consumers",
                                            "coalexitCompensation for Investors Treatment" = "Compensation for Investors",
                                            "coalexitCompensation for Regions" = "Compensation for Regions",
                                            "coalexitCompensation for Workers" = "Compensation for Workers",
                                            "party_2021_2" = "Party Vote 2021",
                                            "left_right_4" = "Left-Right Leaning",
                                            "climate_policyn"= "Climate Policy Support",
                                            "climatehumans" = "Belief in MM CC",
                                            "highintervention" = "Belief in St. Int. Bin.",
                                            "supportintervention" = "Belief in St. Int. Cont.",
                                            "region" = "Region "),
                            estimate  = "{estimate} ({std.error})",
                            statistic = NULL,
                            output="latex")
  
  table.rf2 
  
  regtable.rf2 <- table.rf2 %>% kableExtra::add_header_above(header=tableheader.tr1) 
  
  
  save_kable(regtable.rf2, paste0("Tables/coal_vignette_robust2.weighted.,",file.names.y[i],".tex"))
  
  
  
  
  
}


